function [res,TP] = getTermPremia(model,x_cu)

%% Computing Mean and standard deviation for term premia
% The fitted value for term premium
T            = size(x_cu,2);
TP           = zeros(length(model.TP0),T);
for t=1:T
    if model.pruningOn == 1
        if model.appOrder == 1
            TP(:,t) = model.TP0 + model.TPx*x_cu(1:model.nx,t);
        elseif model.appOrder == 2
            xf_cu                 = x_cu(1:model.nx,t);
            xs_cu                 = zeros(model.nx,1);
            xs_cu(1:model.nx1,1)  = x_cu(model.nx+1:model.nx+model.nx1,1);
            TP(:,t) = model.TP0 + model.TPx*(xf_cu+xs_cu) + model.TPxx*kron3(xf_cu,xf_cu);
        elseif model.appOrder == 3
            xf_cu                 = x_cu(1:model.nx,t);
            xs_cu                 = zeros(model.nx,1);
            xs_cu(1:model.nx1,1)  = x_cu(model.nx+1:model.nx+model.nx1,t);
            xrd_cu                = zeros(model.nx,1);
            xrd_cu(1:model.nx1,1) = x_cu(model.nx+model.nx1+1:model.nx+2*model.nx1,t);
            AA_cu                 = kron3(xf_cu,xf_cu);
            AAA_cu                = kron(xf_cu,AA_cu);
            BB_cu                 = kron3(xf_cu,xs_cu);
            TP(:,t)   = model.TP0 + model.TPx*(xf_cu+xs_cu+xrd_cu) + model.TPxx*(AA_cu+2*BB_cu) ...
                                   + model.TPxxx*AAA_cu;
        end
    else
        AA_cu   = kron3(xhat(1:model.nx,t),xhat(1:model.nx,t));
        TP(:,t) = model.TP0 + model.TPx*xhat(1:model.nx,t) + model.TPxx*AA_cu + model.TPxxx*kron3(AA_cu,xhat(1:model.nx,t));
    end
end

%% The output - multiplying by 40000 (to get tp in annual basis points)
res.MeanTP     = mean(TP*40000,2);
res.stdTP      = std(TP*40000,0,2);
if T < 1000
    res.TP     = TP*40000;
end


end